The distribution of local fluxes in porous media 
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We study the distributions of channel openings, local fluxes, and velocities in a two-dimensional 
random medium of non-overlapping disks. We present theoretical arguments supported by numeri- 
cal data of high precision and find scaling laws as function of the porosity. For the channel openings 
we observe a crossover to a highly correlated regime at small porosities. The distribution of veloc- 
ities through these channels scales with the square of the porosity. The fluxes turn out to be the 
convolution of velocity and channel width corrected by a geometrical factor. Furthermore, while 
the distribution of velocities follows a Gaussian, the fluxes are distributed according to a stretched 
exponential with exponent 1/2. Finally, our scaling analysis allows to express the tortuosity and 
pore shape factors from the Kozeny-Carman equation as direct average properties from microscopic 
quantities related to the geometry as well as the flow through the disordered porous medium. 
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Fluid flow through a porous medium is of importance 
in many practical situations ranging from oil recovery to 
chemical reactors and has been studied experimentally 
and theoretically for a long time QJ@]- Due to disorder, 
porous media display many interesting properties that 
are however difficult to handle even numerically. One 
important feature is the presence of heterogeneities in 
the flux intensities due the varying channel widths. They 
are crucial to understand stagnation, filtering, dispersion 
and tracer diffusion. These are subjects of much practical 
interest in medicine, chemical engineering and geology 
and on which a vast literature is available Q . 

Many stochastic models for disordered porous media 
have been proposed and used to describe the above men- 
tioned effects. One of the most successful is the so-called 
q-model for force distributions in random packings in 
which a scalar fluid is transfered downwards from layer 
to layer. Although the distribution of local flux intensi- 
ties should be the basis for any quantitative evolution of 
these stochastic models, detailed studies of them at the 
pore level are still lacking. 

The traditional approach for the investigation of single- 
phase fluid flow at low Reynolds number in disordered 
porous media is to characterize the system in terms of 
Darcy's law 0, which assumes that a macroscopic 
index, the permeability K, relates the average fluid ve- 
locity V through the pores with the pressure drop AP 
measured across the system, 
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K AP 
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(1) 



where L is the length of the sample in the flow direction 
and n is the viscosity of the fluid. In fact, the permeabil- 
ity reflects the complex interplay between porous struc- 
ture and fluid flow, where local aspects of the pore space 
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morphology and the relevant mechanisms of momentum 
transfer should be a deq uately considered. In previous 
studies 0, S S S S Il0> LUl , computational simulations 
based on detailed models of pore geometry and fluid flow 
have been used to predict permeability coefficients as well 
as to validate semi-empirical correlations obtained from 
real porous materials. 

In this paper we present numerical calculations for a 
fluid flowing through a two-dimensional channel of width 
L y and length L x filled with randomly positioned circu- 
lar obstacles. For instance, this type of model has been 
frequently used to study flow through fibrous filters 01 • 
Here the fluid flows in the x-direction at low but non- 
zero Reynolds number and in the y-direction we impose 
periodic boundary conditions. We consider a particular 
type of random sequential adsorption (RSA) model 0] 
in two dimensions to describe the geometry of the porous 
medium. As shown in Fig. 1, disks of diameter D are 
placed randomly by first choosing from a homogeneous 
distribution between D/2 and L x — D/2 (L y — D /2) the 
random a:-(y-)coordinates of their center. If the disk al- 
located at this position is separated by a distance smaller 
than D/10 or overlaps with an already existing disk, this 
attempt of placing a disk is rejected and a new attempt is 
made. Each successful placing constitutes a decrease in 
the porosity (void fraction) e by ttD 2 / AL x L y . One can 
associate this filling procedure to a temporal evolution 
and identify a successful placing of a disk as one time 
step. By stopping this procedure when a certain value 
of e is achieved, we can produce in this way systems of 
well controlled porosity. We study in particular configu- 
rations with e = 0.6, 0.7, 0.8 and 0.9. 

First we analyze the geometry of our random config- 
urations making a Voronoi construction of the point set 
given by the centers of the disks 0, 0] . We define two 
disks to be neighbors of each other if they are connected 
by a bond of the Voronoi tessellation. These bonds con- 
stitute therefore the openings or pore channels through 
which a fluid can flow when it is pushed through our 
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FIG. 1: Contour plot of the velocity magnitude for a typical 
realization of a pore space with porosity e = 0.7 subjected to 
a low Reynolds number and periodic boundary conditions ap- 
plied in the ^-direction. The fluid is pushed from left to right. 
The colors ranging from blue (dark) to red (light) correspond 
to low and high velocity magnitudes, respectively. The close- 
up shows a typical pore opening of length I across which the 
fluid flows with a line average velocity v. The local flux at 
the pore opening is given by q = vl cos 8, where 9 is the angle 
between v and the vector normal to the line connecting the 
two disks. 




FIG. 2: Distributions of the normalized channel widths I* = 
l/D for different values of porosity e. From left to right, the 
two vertical dashed lines indicate the values of the minimum 
distance between disks I* —0.1 and the size of the disks I* = 
1. The inset shows the data collapse obtained by rescaling 
the distributions with (I) using Eq. 



porous medium, as can be seen in the close-up of Fig. 1. 
We measure the channel widths I as the length of these 
bonds minus the diameter D and plot in Fig. 2 the (nor- 
malized) distributions of the normalized channel widths 
I* = l/D for the four different porosities. Clearly one no- 
tices two distinct regimes: (J) for large widths /* the dis- 
tribution decays seemingly exponentially with I*, and (ii) 
for small I* it has a strong dependence on the porosity, in- 
creasing dramatically at the origin with decreasing poros- 
ity. A closer investigation shows that in Fig. 2 the large 
I* tail decays like a Gaussian for large porosities while 
it is a simple exponential when the porosity is around 
or below 0.7. The crossover between the two regimes is 
visible as a peak which shifts between e = 0.9 and 0.8 
and then stays for smaller porosities at about I* = 1, 
i.e., I — D. These distribution functions can be quali- 
tatively understood in the following way. For very large 
porosities, i.e., very dilute systems, the distance between 
the particles is essentially uncorrelated due to excluded 
volume and is therefore Gaussian distributed around a 
mean value (I). If for simplicity one imagines particles 
being on a regular triangular lattice as an idealized con- 
figuration in two dimensions, the following expression is 
obtained: 

(|)=g(./II 7r -1). (2) 
W V V 2V3(1 -e) ' 1 ' 

The filling process will strongly feel the clogging due 
to excluded volume when one disk just fits into the 
hole between three disks. This situation occurs when 
(I) = D(y/3 - 1). Inserting this into Eq. 10) gives a 
crossover porosity of e = 1 — 7r/6v3 « 0.7 which agrees 
with our simulation (see Fig. 2). Interestingly, a related 



property, namely the correlation function, does not seem 
to show such a crossover 0, |^. The inset of Fig. 2 
shows that, for sufficiently large values of I, all distribu- 
tions P(l) collapse to a single curve when rescaled by the 
corresponding value of (I) calculated from Eq. J2J. As 
shown in the inset of Fig. 3, the variation of the aver- 
age value (I*) with porosity follows very closely Eq. J2J. 
Only the prefactor is different from unity (w 1.2) due to 
the presence of disorder. This result indicates that our 
simple description based on a diluted system of particles 
placed on a regular lattice provides a good approximation 
for the geometry of the disordered porous medium. 

The fluid mechanics in the porous space is based on the 
assumption that a Newtonian and incompressible fluid 
flows under steady-state conditions. The Navier-Stokes 
and continuity equations for this case reduce to 

p u-Vu = -Vp + pV 2 u , (3) 



V • u = , (4) 

where u and p are the local velocity and pressure fields, 
respectively, and p is the density of the fluid. No- 
slip boundary conditions are applied along the entire 
solid-fluid interface, whereas a uniform velocity profile, 
u x (0, y) — V and u y (0,y) = 0, is imposed at the in- 
let of the channel. For simplicity, we restrict our study 
to the case where the Reynolds number, defined here as 
Re = pVLy/p, is sufficiently low (Re < 1) to ensure a 
laminar viscous regime for fluid flow. We use FLUENT 
|l8j . a computational fluid dynamic solver, to obtain the 
numerical solution of Eqs. © and on a triangulated 
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FIG. 3: Double-logarithmic plot of the distributions of the 
local normalized velocity magnitudes v* , i.e., v/V , multiplied 
by e 2 as explained in the text. The solid line is a Gaussian 
fit. The inset shows the dependence of (l*) and (v*) on the 
porosity e. The solid lines are the best fits to the data, cor- 
responding to (P) = a(^/b(l — e) — 1), with a — 1.22 and 
b = tt/2^3 (see Eq. and (v*) = 0.71e -2 . 



grid of up to hundred thousand points adapted to the 
geometry of the porous medium. 

Simulations have been performed by averaging over 10 
different pore space realizations generated for each value 
of porosity. The contour plot in Fig. 1 of the local ve- 
locity magnitude for a typical realization of the porous 
medium with porosity e = 0.7 clearly reveals that the 
transport of momentum through the complex geometry 
generates preferential channels [Tl| . Once the numerical 
solution for the velocity and pressure fields in each cell 
of the numerical grid is obtained, we compute the fluid 
velocity magnitudes v associated to each channel. This 
value is the magnitude of the line average velocity vector 
v calculated as the average over the local velocity vectors 
u along the corresponding channel width I. 

In Fig. 3 we show the data collapse of all distribu- 
tions of normalized velocity magnitudes P(v*), where 
v* = v/V, obtained by rescaling the variable v* with the 
corresponding value of e~ 2 . It is also interesting to note 
that these rescaled distributions follow a typical Gaus- 
sian behavior except for very small v*e 2 , as indicated by 
the solid line in Fig. 3. In the inset of Fig. 3 we also 
show that the average interstitial velocity indeed scales 
with the porosity as (v) ~ e~ 2 , confirming the rescaling 
procedure adopted to obtain the collapse of the distribu- 
tions P(v*) in the main plot of Fig. 3. Plotting for each 
channel v against I gives a cloud of points which for all 
considered values of e results in a rather unexpected least 
square fit relation of the type v ~ \fl. 

We now analyse the distribution of fluxes throughout 
the porous medium. Each local flux q crossing its corre- 
sponding pore opening I is given by q = vl cos 9, where 9 



is the angle between v and the vector normal to the cross 
section of the channel (see Fig. 1). In Fig. 4 we show that 
the distributions of normalized fluxes cf> = q/qt, where 
qt = VL y is the total flux, have a stretched exponential 
form, 



P(4>) ~ exp(-y/j/fa) 



(5) 



with c/)q w 0.005 being a characteristic value. This sim- 
ple form of Eq. © is quite unexpected considering the 
rather complex dependence of P(l) on e. Moreover, all 
flux distributions P((f>) collapse on top of each other when 
rescaled by the corresponding value of e 2 . This col- 

lapse for distinct porous media results from the fact that 
mass conservation is imposed at the microscopic level of 
the geometrical model adopted here, which is microscop- 
ically disordered, but at a larger scale is macroscopically 
homogeneous Q. As also shown in Fig. 4, it is possible 
to reconstruct the distribution of fluxes using a convolu- 
tion of the distribution of velocities v and the distribution 
of oriented channel widths, namely lcos9. Indeed, if we 
calculate the integral, 



P{4>) = J J P(v)P(l cos 9)6(<f>-vl cos 6)dvd{l cos 6) , (6) 

we find that the original distribution P(4>) is approxi- 
mately retrieved, as can also be seen in Fig. 4 (solid line). 

Finally, the inset of Fig. 4 shows that the permeability 
of the two-dimensional porous media closely follows the 
semi-empirical Kozeny-Carman equation 0, 



K 



1 



(7) 



where K$ = ft, 2 / 12 is a reference value taken as the per- 
meability of an empty channel between two walls sepa- 
rated by a distance h. The proportionality constant k is 
given by the following expression: 



1 

ret 



(8) 



where r = (L e /L) 2 is the hydraulic tortuosity of the 
porous medium, a corresponds to the pore shape factor, 
and L e is an effective flow length [l| . If we now make use 
of the Dupuit-Forchheimer assumption , 
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(9) 



we are led to the conclusion that the tortuosity of our 
porous medium should also scale as t ~ e~ 2 . Considering 
the validity of the Kozeny-Carman equation J7J and the 
definition of the constant k from Eq. JHJl, we obtain as 
a consequence that the shape factor should behave as 
a ~ e 2 . 
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FIG. 4: Log-log plot of the distributions of the normalized 
local fluxes (f> — q/qt for different porosities e. The (red) 
dashed line is a fit of the form exp(— \J (f>/(f>o), where cpo ~ 
0.005. The full line stems from the convolution as discussed 
in the text. In the inset we see a double-logarithmic plot of the 
global flux and the straight line verifies the Kozeny-Carman 
equation. 



Summarizing we have found that although the distribu- 



tion of channel widths in a porous medium made by a 
two-dimensional RSA process is rather complex and ex- 
hibits a crossover at I ~ D, the distribution of fluxes 
through these channels shows an astonishingly simple be- 
havior, namely a square-root stretched exponential dis- 
tribution that scales in a simple way with the porosity. 
The velocity magnitudes follow a Gaussian distribution 
truncated at small velocities which scales with the square 
of the porosity. The distribution of fluxes can be recon- 
structed as a convolution of the velocity with the channel 
widths distributions corrected by the velocity orientation 
factor cos 9. We propose simple scaling laws for the local 
fluxes that deepen the understanding of the intrinsic con- 
nection between geometrical and flow properties of the 
random porous medium. Furthermore, we show that our 
results can be macroscopically described in terms of the 
Kozeny-Carman relation. Future tasks consist in gener- 
alizing these studies to higher Reynolds numbers, three 
dimensional model of porous media and other types of 
disorder. Other important challenges are to investigate 
transient flow and tracer dynamics. 
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